上一篇介绍了显式欧拉法、隐式欧拉法、两步欧拉法和改进欧拉法求解常微分方程初值问题；其中显式欧拉法和隐式欧拉法是一阶算法精度，截断误差为 O ( h 2 ) O\left( {{h^2}} \right) O(h2)；两步欧拉法和改进欧拉法是二阶算法精度，截断误差为 O ( h 3 ) O\left( {{h^3}} \right) O(h3)；欧拉法的精度有限、需要求解步长 h h h很小。本篇介绍求解精度更高的四阶龙格库塔法(Runge-Kutta)，其截断误差为 O ( h 5 ) O\left( {{h^5}} \right) O(h5)。

对于一阶微分方程初值问题：

{ y ˙ = f ( t , y ) y ( t 0 ) = y 0 \left\{ \begin{array}{l} {\bf{\dot y}} = f\left( {t,{\bf{y}}} \right) \\ {\bf{y}}\left( {{t_0}} \right) = {{\bf{y}}_0} \\ \end{array} \right. {y˙​=f(t,y)y(t0​)=y0​​

式中， t 0 {t_0} t0​为初始时间（已知常数）， y 0 {{\bf{y}}_0} y0​为初始状态（已知向量）， f ( t , y ) f\left( {t,{\bf{y}}} \right) f(t,y)为关于时间 t {t} t和状态 y {{\bf{y}}} y的函数（已知函数）。

四阶龙格库塔法(Runge-Kutta)求解算法为：

k 1 = f ( t ( k ) , y ( k ) ) {{k}_{1}}=f\left( t\left( k \right),\mathbf{y}\left( k \right) \right) k1​=f(t(k),y(k)) k 2 = f ( t ( k ) + h 2 , y ( k ) + h 2 k 1 ) {{k}_{2}}=f\left( t\left( k \right)+\frac{h}{2},\mathbf{y}\left( k \right)+\frac{h}{2}{{k}_{1}} \right) k2​=f(t(k)+2h​,y(k)+2h​k1​) k 3 = f ( t ( k ) + h 2 , y ( k ) + h 2 k 2 ) {{k}_{3}}=f\left( t\left( k \right)+\frac{h}{2},\mathbf{y}\left( k \right)+\frac{h}{2}{{k}_{2}} \right) k3​=f(t(k)+2h​,y(k)+2h​k2​) k 4 = f ( t ( k ) + h , y ( k ) + h k 3 ) {{k}_{4}}=f\left( t\left( k \right)+h,\mathbf{y}\left( k \right)+h{{k}_{3}} \right) k4​=f(t(k)+h,y(k)+hk3​) y ( k + 1 ) = y ( k ) + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) \mathbf{y}\left( k+1 \right)=\mathbf{y}\left( k \right)+\frac{h}{6}\left( {{k}_{1}}+2{{k}_{2}}+2{{k}_{3}}+{{k}_{4}} \right) y(k+1)=y(k)+6h​(k1​+2k2​+2k3​+k4​) y ( 0 ) = y 0 \mathbf{y}\left( 0 \right)={{\mathbf{y}}_{0}} y(0)=y0​

上式是关于 y ( k ) {\bf{y}}\left( k \right) y(k)向 y ( k + 1 ) {\bf{y}}\left( k+1 \right) y(k+1)的递推形式，可以根据初始条件按照递推关系依次求出 y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) ⋯   , y ( N ) ⋯ {\bf{y}}\left( 1 \right),{\bf{y}}\left( 2 \right),{\bf{y}}\left( 3 \right),{\bf{y}}\left( 4 \right) \cdots ,{\bf{y}}\left( N \right) \cdots y(1),y(2),y(3),y(4)⋯,y(N)⋯，此离散序列即为微分方程的数值解。

微分方程的数值解法，本质是使用数值积分来实现 y ˙ {\bf{\dot y}} y˙​向 y {\bf{y}} y的转换。四阶龙格库塔法通过对微分的四步分段逼近，在一个求解步长内能够逼近复杂的曲线，因此能够取得较高的计算精度，其截断误差为 O ( h 5 ) O\left( {{h^5}} \right) O(h5)。

作者使用Matlab开发了四阶龙格库塔法求解常微分方程的程序，能够方便快捷的求解一阶常微分方程的初值问题。

下面结合实例进行演示和分析。

求解一阶常微分方程（式中向量 x {\bf{x}} x等价于前文中的向量 y {\bf{y}} y）： x ˙ = f ( t , x ) = [ x ( 2 ) ∗ x ( 3 ) − x ( 1 ) ∗ x ( 3 ) − 0.51 ∗ x ( 1 ) ∗ x ( 2 ) ] \mathbf{\dot{x}}=f\left( t,\mathbf{x} \right)=\left[ \begin{matrix} \mathbf{x}(2)*\mathbf{x}(3) \\ -\mathbf{x}(1)*\mathbf{x}(3) \\ -0.51*\mathbf{x}(1)*\mathbf{x}(2) \\ \end{matrix} \right] x˙=f(t,x)=⎣⎡​x(2)∗x(3)−x(1)∗x(3)−0.51∗x(1)∗x(2)​⎦⎤​ x ( 0 ) = [ 0 1 1 ] \mathbf{x}\left( 0 \right)=\left[ \begin{matrix} 0 \\ 1 \\ 1 \\ \end{matrix} \right] x(0)=⎣⎡​011​⎦⎤​

不同时间步长 h h h时的数值计算结果：

步长 h = 0.6 s h=0.6s h=0.6s

步长 h = 0.9 s h=0.9s h=0.9s

从计算结果可以看出，四阶龙格库塔法(Runge-Kutta)即使在步长很大时，也能保持较高的求解精度，求解精度与Matlab自带的ode45函数相当，相对于改进欧拉算法求解精度有明显提高。 自己编程实现四阶龙格库塔法(Runge-Kutta)，相对于直接调用ode45等Matlab自带的龙格库塔法的最大优势在于：可以将求解程序和模型描述文件融合起来，解决各类参数时变、条件判断、多模型切换等问题。后续补充实际案例。

有Matlab/Simulink方面的技术问题，欢迎发送邮件至[email protected]讨论。 更多Matlab/Simulink原创资料，欢迎关注微信公众号：Matlab Fans

源程序下载: 1. csdn资源: 四阶龙格库塔法(Runge-Kutta)求解常微分方程的Matlab程序及案例 2. 扫码关注微信公众号Matlab Fans，回复BK09获取百度网盘下载链接。